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Resumen 


Aragón-Hernández, J. L., : Bladé, E. (mayo-junio, 2017). 
Modelación numérica de flujo mixto en conductos cerrados 
con esquemas en volúmenes finitos. Tecnología y Ciencias del 
Agua, 8(3), 127-142. 


En este trabajo se presenta un modelo numérico para la 
simulación de flujo mixto (flujo en lámina libre y flujo en 
presión) en conductos cerrados a través de las ecuaciones de 
Saint Venant en una dimensión para flujo en lámina libre y 
con el método de la ranura de Preissmann para considerar el 
flujo en presión. Para la resolución de las ecuaciones se emplea 
el método de los volúmenes finitos con un esquema de alta 
resolución. El esquema utilizado es el método de Godunov 
con el Riemann solver de Roe de primer orden de precisión, 
más unas correcciones de segundo orden para obtener el 
esquema de alta resolución WAF-TVD. El flujo mixto es un 
fenómeno bastante común en colectores pluviales, túneles, 
tuberías de obras de toma de instalaciones hidroeléctricas, 
llenado / vaciado de tuberías, colectores de almacenamiento, 
etcétera. La entrada en carga en los conductos se puede 
generar desde el extremo aguas abajo, desde el extremo 
aguas arriba y por ambos extremos simultáneamente, siendo 
la primera la más común y, por lo tanto, la evaluada en este 
trabajo. Para demostrar la actuación del modelo, éste se 
aplica a un caso de referencia y a dos ensayos de laboratorio 
existentes en la literatura técnica. Los resultados obtenidos 
muestran que el modelo numérico es capaz de reproducir 
los experimentos con buena precisión. Con ello se demuestra 
que el modelo es idóneo para simular flujo en lámina libre en 
régimen lento, rápido y transcrítico (de lento a rápido y de 
rápido a lento), y flujo mixto. 


Palabras clave: flujo en lámina libre, flujo en presión, flujo 
mixto, ranura de Preissmann, volúmenes finitos y colectores 
pluviales. 


Abstract 


Aragón-Hernández, ]. L., € Bladé, E. (May-June, 2017). Numerical 
modeling of mixed flow in closed conduits with finite volume 
schemes. Water Technology and Sciences (in Spanish), 8(3), 
127-142. 


This paper presents a numerical model for simulating mixed flow 
(free surface and pressurized flow) in closed conduits through the 
one-dimensional Saint Venant equations for free surface flow and the 
method Preissmann slot to consider the pressurized flow. To solve the 
equations the finite volume method with a high resolution scheme is 
used. The used scheme is the Godunov method with the Roe solver 
Riemann with precision of first order, plus a second-order corrections 
to obtain the high resolution WAF-TVD scheme. The mixed flow 
is a common phenomenon in storm sewers, tunnels, pipes of 
hydroelectric installations, filling/draining pipes, storage collectors, 
etc. The pressurized flow in close conduits can be generated from 
downstream end, upstream end and both ends simultaneously; the 
first one is the most common, therefore will be evaluated in this 
study. To demonstrate the performance of the model, the numerical 
model is applied to a reference test and a two laboratory tests existing 
in the technical literature. The results show that the numerical model 
is able to reproduce the experiments with good accuracy. This shows 
that the model is suitable to simulate subcritical, supercritical and 
transcritical (from subcritical to supercritical and from supercritical 
to subcritical flow) free surface and mixed flow. 


Keywords: Free surface flow, pressurized flow, mixed flow, 
Preissmann slot, finite volumes and storm sewers. 
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Introducción 


La ocurrencia simultánea de flujo en lámina 
libre y flujo en presión en un conducto cerrado 
se conoce comúnmente como flujo mixto. Este 
fenómeno se puede encontrar de manera fre- 
cuente en colectores pluviales, túneles, tuberías 
de obras de toma de instalaciones hidroeléctri- 
cas, llenado /vaciado de tuberías, colectores 
de almacenamiento, etcétera. El proceso de 
transición de flujo en lámina libre a flujo en 
presión puede ser inducido, entre otras causas, 
por la variación de caudales de entrada, bloqueo 
o reducción de la capacidad de la tubería, su- 
mersión de la salida de la tubería, operación de 
elementos de control (vertedores, compuertas, 
etc.), fallo de estaciones de bombeo y presencia 
de pozos de caída (Bourdarias éz Gerbi, 2007; Li 
€ McCorquodale, 1999). 

En el ámbito urbano, el rol tradicional de una 
red de alcantarillado ha sido principalmente la 
conducción de aguas pluviales y aguas residua- 
les. La gran mayoría de los colectores pluviales 
de una red de drenaje se diseñan para operar en 
condiciones normales con flujo en lámina libre y 
régimen lento, considerando caudales de diseño 
para eventos de precipitación bajos. Cuando la 
capacidad de uno o más colectores es excedida 
por un evento de precipitación, lo que no es 
raro, O por cambios bruscos en las condiciones 
de contorno, el tipo de flujo puede cambiar de 
lámina libre a flujo en presión. 

Durante la transición de flujo en lámina 
libre a flujo en presión pueden ocurrir fuertes 
transitorios de presión y provocar daños a los 
colectores, levantamiento de tapas de los pozos 
de visita, expulsión de agua e inundación de la 
superficie urbana (Li éÉ McCorquodale, 1999). 
Los frentes de onda generados en la transición 
de flujo en lámina libre a flujo en presión pue- 
den ser de dos tipos: gradual, que consiste de un 
frente de onda en lámina libre con condiciones 
de flujo en lámina libre en ambos lados del fren- 
te, y brusco, que es un frente de onda a conducto 
lleno con condiciones de flujo en lámina libre 
adelante y de flujo en presión atrás del frente 
(Wright, Vasconcelos, Creech, $1 Lewis, 2008). 


Para modelar el flujo mixto existen diferentes 
aproximaciones, como el método de la columna 
rígida (Li é£ McCorquodale, 1999) y las incluidas 
dentro de los modelos dinámicos. Los modelos 
dinámicos se pueden clasificar en aquellos que 
utilizan uno o dos sistemas de ecuaciones. Los 
métodos que utilizan un sistema de ecuacio- 
nes requieren de una aproximación adicional 
para considerar flujo en presión. Entre estas 
aproximaciones se encuentran el método TPA 
(Two-Component Pressure Approach) (Vasconcelos, 
Wright, £: Roe, 2006) y el método de la ranura 
de Preissmann, uno de los más conocidos y 
utilizados. 

Preissmann (1961) fue quien reconoció la 
similitud de la forma entre canales abiertos y 
conductos cerrados; además, propuso que una 
ranura hipotética se agregara a la clave de los 
conductos cerrados, tal que tanto el flujo en 
lámina libre como el flujo en presión pudieran 
ser analizados mediante las ecuaciones de Saint 
Venant en una dimensión (Li £ McCorquodale, 
1999). El método de la ranura de Preissmann 
fue implementado por primera vez por Cunge 
y Wegner (1964), al estudiar flujo en presión, 
posteriormente varios autores (Capart, Sillen, éz 
Zech, 1997; García-Navarro, Priestley, €: Alcru- 
do, 1994; Ji, 1998; Kerger, Archambeau, Erpicum, 
Dewals, éz Pirotton, 2011; León, Ghidaoui, Sch- 
midt, 6: García, 2009; Song, Cardle, € Leung, 
1983; Trajkovic, Ivetic, Calomino, € D'Ippolito, 
1999) han desarrollado modelos numéricos ba- 
sados en este método para simular flujo mixto. 

Cuando se usa esta aproximación pueden 
aparecer cuatro limitaciones principales aso- 
ciadas con la aproximación misma y con el 
esquema numérico empleado para resolver las 
ecuaciones del flujo. La primera es que asume 
flujo a superficie libre en todo el sistema; por 
consiguiente, si la carga piezométrica cae por 
debajo de la clave del colector se desarrolla 
flujo en lámina libre y por tanto es incapaz de 
simular flujo en presión subatmosférica (León 
et al., 2009, León, Ghidaoui, Schmidt, éz García, 
2010; Li £e McCorquodale, 1999; Politano, Odga- 
ard, éz Klecan, 2007; Vasconcelos et al., 2006); la 
segunda, inestabilidades asociadas con el pobre 
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desempeño del esquema numérico cuando el 
flujo cambia rápidamente de flujo en lámina li- 
bre a flujo en presión, la tercera, inconsistencias 
en el balance de masa y momento asociados con 
el ancho de la ranura; la cuarta, inexactitudes 
en la propagación de transitorios hidráulicos 
relacionados con el ancho de la ranura (León 
et al., 2009). 

Para superar la primera limitación, reciente- 
mente Kerger et al. (2011) ampliaron el concepto 
de la ranura de Preissmamn, desarrollando la 
ranura de Preissmann negativa, que consiste 
en extender la ranura por debajo de la clave 
de la tubería para modelar flujos con presiones 
subatmosféricas o negativas, y aire atrapado, 
acoplando un modelo matemático de equilibrio 
homogéneo. Con respecto a la segunda, León et 
al. (2009) propusieron una modificación al mé- 
todo, introduciendo una transición geométrica 
gradual entre la tubería y la ranura, desde 0.95 
hasta 1.5 veces el diámetro de la tubería. Esta 
modificación permite que exista continuidad 
en la velocidad de la onda en lámina libre y 
en presión, en función del ancho de la ranura; 
mientras que Barnett (2010) sugiere una transi- 
ción geométrica exponencial del tipo w = Te*; 
w es el ancho de la ranura al nivel z, medida 
hacia abajo desde el punto donde w = T, y T, es 
el ancho de la ranura de Preissmann en función 
de la celeridad de la onda de presión a. La forma 
de la ranura depende esencialmente del pará- 
metro k; valores más pequeños corresponden 
a transiciones más elongadas. Con este tipo de 
transiciones existe un aumento del área de la 
sección tanto en el colector como en la ranura, 
que conlleva un almacenamiento extra no real 
en dicha transición, por lo que pueden aumen- 
tar los problemas relacionados con la tercera y 
cuarta limitación. 

El método de la ranura de Preissmann ha 
sido usado para la modelación de transición 
gradual (suave) de flujo mixto con éxito, y 
transición brusca (rápida) con poco éxito o 
sin éxito. El poco éxito en la modelación es 
principalmente consecuencia de los problemas 
de inestabilidad si la onda de presurización es 
demasiado brusca, lo cual puede causar que 


la simulación aborte incluso para anchos de la 
ranura 0.1 veces el diámetro (León et al., 2009). 

Por otra parte, para la solución de sistemas 
de ecuaciones hiperbólicas, como las ecuaciones 
de Saint Venant, actualmente existen diferentes 
esquemas, entre ellos los métodos directos 
(Shock Capturing Methods), que son muy atrac- 
tivos, pues permiten describir o capturar dis- 
continuidades del flujo (p. ej., frentes de onda, 
resaltos hidráulicos) de forma automática (Toro, 
1999). Estos métodos inicialmente fueron desa- 
rrollados en el campo de la dinámica de gases 
y después extendidos a la hidráulica de ríos y 
canales (Bladé, Gómez-Valentín, Sánchez-Juny, 
éz Dolz, 2008; Toro éz García-Navarro, 2008). En 
los últimos años, estos esquemas han sido uti- 
lizados para el cálculo de flujo mixto (Aragón- 
Hernández, Concha, Bladé, $2 Gómez, 2009; 
Capart et al., 1997; León et al., 2009; Sanders éz 
Bradford, 2011; Vasconcelos et al., 2006). 

El objetivo de este trabajo es presentar el 
desarrollo de un modelo numérico estable y 
robusto, incluso para transiciones bruscas, para 
reproducir flujo mixto en conductos cerrados. 
El modelo emplea las ecuaciones completas de 
Saint Venant en una dimensión para el flujo 
en lámina libre y el método de la ranura de 
Preissmann para considerar flujo en presión. 
La resolución numérica de dichas ecuaciones se 
realiza con el método de los volúmenes finitos 
con un esquema conservativo y descentrado de 
primer orden más unas correcciones de segundo 
orden, para finalmente obtener un esquema con 
una precisión de alta resolución. En específico, 
este esquema numérico ampliamente utilizado 
en ríos será evaluado por primera vez al mo- 
delar flujo mixto en conductos cerrados y de 
entrada permitirá superar la segunda limitación 
mencionada. 

Este modelo forma parte de una investiga- 
ción más amplia relacionada con la modelación 
numérica de la dinámica del agua en zonas 
urbanas. El modelo se integra dentro de un 
modelo más general que se utiliza para simular 
el comportamiento hidráulico del flujo de agua 
en el sistema de drenaje superficial o sistema de 
drenaje mayor (canales, tejados, calles, aceras) 
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y en el sistema de drenaje subsuperficial o 
sistema de drenaje menor (colectores, unión de 
colectores), integrado completamente a través 
del intercambio de flujo en los puntos físicos 
correspondientes (rejas, pozos de visita) en am- 
bas direcciones, conocido como drenaje urbano 
dual (Aragón-Hernández, 2013). 

Para validar el modelo numérico desarro- 
llado, se emplea un caso de referencia y dos 
ensayos de laboratorio existentes en la literatura 
técnica. En dichos casos de validación se incluye 
la mayoría de los tipos de flujo que de forma 
general se pueden presentar en conductos cerra- 
dos y, sobre todo, en colectores pluviales, como 
flujo en lámina libre en régimen lento, rápido y 
transcrítico, y flujo mixto. 

Así, el trabajo se estructura en cinco aparta- 
dos. En el primero se describen las ecuaciones 
para flujo en lámina libre; en el segundo, la 
aproximación para considerar flujo en presión; 
en el tercero se plantea el esquema numérico; 
en el cuarto se hace la validación numérica, y 
finalmente, en el quinto se presentan las con- 
clusiones. 


Flujo en lámina libre 


Para modelar flujo en lámina libre se utilizan las 
ecuaciones completas de Saint Venant en una 
dimensión. La deducción de estas ecuaciones 
se realiza de manera directa aplicando las leyes 
de conservación de la masa y cantidad de movi- 
miento a un volumen de control, considerando 
una sección arbitraria y canal no prismático, 
resultando un sistema de ecuaciones en forma 
conservativa (Chaudhry, 1993): 


PA (1) 
dt  0x 

9Q,Ia[(Q 

22 ]as,+5(8,-5) (2) 


Utilizando notación vectorial resulta: 


Lp. (3) 
ot Ox 


Donde el vector de variables conservativas 
U, el vector de flujo F y el término fuente H 
responden a las expresiones: 


Q 0 


F= 2 H= 
E+8l, 81, +8gA(S,-S, 


[e 
Il 


(4) 
) 


En las expresiones anteriores, A es el área 
de la sección mojada; Q, el caudal circulante; 
g, la aceleración de la fuerza de gravedad, x, 
la coordenada longitudinal en dirección del 
flujo; t, el tiempo; S,, la pendiente del terreno; 
S ñ la pendiente de fricción o motriz; [,, la fuerza 
debida a la presión del agua en una sección o 
momento de primer orden de la sección respecto 
a la superficie libre del agua, e I,es la contribu- 
ción de las fuerzas de presión del contorno; las 
tres últimas se expresan, respectivamente, como: 


L,= f, (M=m)b (x,mjan (6) 
E 
L,= $, (1-0) Ea (7) 


Donde n es el coeficiente de fricción de Man- 
ning; R,, el radio hidráulico; h, el calado; v, la 
velocidad del flujo; b, el ancho de la superficie 
libre del agua, y n es la profundidad del centro 
de gravedad de la sección. Para canales prismá- 
ticos, aunque tengan una sección cualquiera, el 
término l, es idénticamente igual a cero. Tal es 
el caso de la mayoría de los conductos cerrados. 


Flujo mixto 


Como se definió antes, el flujo mixto es la 
ocurrencia al mismo tiempo de flujo en lámina 
libre y flujo en presión en un conducto cerrado. 
Cuando se usa un modelo dinámico de un solo 
sistema de ecuaciones es necesaria una apro- 
ximación adicional; en este trabajo se utiliza el 
método de la ranura de Preissmann. El método 
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consiste en implementar una ranura hipotética 
de ancho T, en la clave de los conductos cerrados 
(figura 1), proporcionando la presión hidrostá- 
tica extra, mientras que se permite el uso de las 
ecuaciones de flujo en lámina libre para ambos 
tipos de flujo (flujo en lámina libre y flujo en 
presión) (Capart, Bogaerts, $: Kevers-Leclercq, 
1999; García-Navarro et al., 1994; León et al., 
2009; Politano et al., 2007; Preissmann, 1961). 

En esencia, esta aproximación explota la 
similitud de la ecuación de onda del flujo en 
lámina libre y flujo en presión. Así, a partir de la 
ecuación (8), que define la celeridad de la onda 
de gravedad, y para evitar caer en el tercer y 
cuarto problema asociados con el ancho de la 
ranura, ésta se calcula bajo la consideración de 
que la celeridad de la onda de presión a es igual 
a la celeridad de la onda de gravedad c y el área 
a tubo lleno A... 


8% 
0 El (8) 


De esta forma, el ancho de la ranura de 
Preissmann T, queda definido como: 


T; > e EE (9) 


Cuando en un colector ocurre flujo en pre- 
sión, el calado calculado es la carga de presión. 
Cuando se usa un ancho de ranura pequeño 
para poder representar grandes niveles de agua 
o una velocidad de onda grande, el método 
es numéricamente inestable (Politano et al., 
2007; Yen, 1978); tales inestabilidades pueden 
ser removidas haciendo ranuras más anchas. 
Sin embargo, cambiar el ancho de la ranura 
destruye la equivalencia con la celeridad de la 
onda mencionada anteriormente. El ancho de la 


Flujo en lámina libre ——> 


Flujo en presión 


Figura 1. Flujo en lámina libre y flujo en presión con el método de la ranura de Preissmamn. 
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ranura es un parámetro significativo y una elec- 
ción óptima del mismo determinará la calidad 
de los resultados. Es importante destacar que al 
implementar este método, la ranura no debería 
aumentar el área de la sección transversal ni 
el radio hidráulico del conducto en presión 
(Chaudhry, 1987), pero en la práctica se permite 
utilizar ranuras con un ancho de hasta 0.1D, 
donde D es el diámetro de la tubería (Trajkovic 
et al., 1999). 

Debido a que se utiliza un único sistema 
de ecuaciones, no hay necesidad de rastrear 
la propagación del movimiento de la interfase 
entre el flujo en lámina libre y flujo en presión 
(Li £ McCorquodale, 1999; Vasconcelos et al., 
2006), el cual puede ser capturado como ondas 
de choque o frentes de onda con los métodos 
directos (Toro, 1999). 


Esquema numérico 


Para la resolución de las ecuaciones de Saint 
Venant en una dimensión se utiliza el método de 
los volúmenes finitos. El método de volúmenes 
finitos se basa en las ecuaciones de gobierno 
escritas en forma integral sobre una celda o 
volumen finito. El dominio de estudio o malla 
de cálculo se discretiza en un dominio espacial x 
en celdas de longitud Ax y un dominio temporal 
t, en intervalos de tiempo Af. Las variables hi- 
dráulicas utilizadas (4, Q) representan el valor 
medio de las variables dependientes en cada 
celda í, centradas en el nodo x, y se extienden 
de ¡- Y a 1 +7. De esta forma, el esquema nu- 
mérico para la solución de la ecuación (3) viene 
dado por: 


n+1 n At sd di At * 
u > u; an Ea) +A (10) 
Donde Uy U”* son los valores medios del 
vector U en la celda ¡ en el instante de tiempo t" y 


t"1, respectivamente; F' es el flujo numérico 


en las interceldas, y H', ein valor representa- 
tivo del promedio del término fuente H en la 
celda ¡ en el instante de tiempo t"y t"". F” ,,, y 
H' en general dependen de las celdas contiguas 


a la celda ¡ en los instantes de tiempo ft" y t”*!, 


El método de Godunov es un esquema con- 
servativo que utiliza una discretización descen- 
trada en volúmenes finitos, cuya particularidad 


es que el flujo numérico F' entre dos celdas i, 


1:1/2 
e 1+ 1 se obtiene de la clan de un problema 
de Riemann local entre dos estados constantes 
U; y U””. El problema de Riemamn tiene una 
estructura compleja y encontrar su solución es 
costoso, por lo cual diversos autores (Harten, 
Lax, €£ Van Leer, 1983; Osher €: Solomon, 1982; 
Roe, 1981) desarrollaron métodos para encon- 
trar una respuesta aproximada del problema de 
Riemanmn, conocidos como aproximate Riemann 
solvers (Toro, 1999). Uno de ellos es el Riemann 
solver de Roe, empleado en este trabajo, el cual 
lleva a un esquema conservativo denominado 
método de Godunov con el Riemann solver 
de Roe; este esquema es de primer orden de 
precisión. 

En zonas cercanas a fuertes gradientes O 
donde se presentan discontinuidades de la 
solución, los esquemas de primer orden son 
disipativos (amplitud de la onda menor a la de 
la solución exacta), mientras que los esquemas 
de segundo orden son dispersivos (producen 
oscilaciones en la solución). Por esta razón, se 
necesitan esquemas de segundo orden o mayor, 
que no presenten dispersión en las discontinui- 
dades, como los esquemas TVD (Total Variation 
Dimishing). Un esquema así se puede obtener 
a partir de un esquema de primer orden más 
un esquema de segundo orden. Para evitar la 
dispersión del esquema de segundo orden se 
debe limitar el flujo numérico del mismo; en este 
trabajo se utiliza el esquema WAF (Weight Avera- 
ged Flux). Así, el esquema numérico WAF-TVD 
resultante se puede escribir como un esquema 
de primer orden más unas correcciones de se- 
gundo orden, para conseguir que la precisión 
sea de alta resolución (Bladé, 2005): 


+[ Duo (poes Bvisteno( Jl) (11) 


ja 
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Donde: 
v=hk—= (12) 


Los coeficientes constantes o las llamadas 
fuerzas de cada onda «:, ,, los valores y vectores 
propios 4, , y £, , de la matriz jacobiana definida 


- gE . ] 
como J = e se escriben, respectivamente, como: 


Ñ ¿EFÚUJAA=AO 
om Ja (13) 
hy, =Úx0 (14) 
loa (15) 


La velocidad í, los incrementos de área 
hidráulica AA y del caudal circulante AQ se 
definen como: 


Q, + Qi 
ú= JA, Aja (16) 
(157+/a,,) 
AA=A, A, (17) 
AQ=Q,.4-Q, (18) 


De la descomposición de los vectores propios 
resulta Y, ,, definida como: 


A 1 
Yi =>584(14) (19) 
Y ==Y 
La corrección de entropía de Harten e 


Hyman q, según Toro (1999), es: 


P, 


á 28€ 
Sl 1 


E, si 0 


Donde: 


e, im/2 A máx [0,7 ee Das Jal Aia = 3] (21) 


Algunas formulaciones utilizadas para la 
celeridad promedio í no reproducen el salto de 
las fuerzas de presión, por lo que en este trabajo 
se utiliza una expresión que mantiene el signi- 
ficado físico de la celeridad, entendido como la 
variación de las fuerzas de presión respecto a 
la variación del área de flujo en cada sección 
transversal (Bladé et al., 2008): 


1 -Ij-—(AL|. 
mE e) ( aba a: Apt 
¡+1 i (22) 
si A, =Aja 


El área hidráulica Á y el ancho de la superfi- 
cie libre del agua B son: 


Á = 2AA, 41 (23) 
A, +A 
BE B, Ss (24) 


La variación de las fuerzas de presión 
(Az, | las para un área constante A, ,, es: 


(AD. L4).., = (Bla). -(2,/4), (25) 


Donde A) y aL, son las fuerzas de 
presión en el volumen finito ¡ e ¿+ 1 para un 
valor del área A. 

Una de las funciones de limitación Y; más 
utilizadas y también la empleada en este trabajo 


es la función la de Minmod, definida como: 


y, = máx [0,mín (1,1)] (26) 


I55N 0187-8336 + 


Tecnología y Ciencias del Agua, vol. VU, núm. 3, mayo-junio de 2017, pp. 127-142 


Tecnología y Ciencias del Agua, vol. VHI, núm. 3, mayo-junio de 2017, pp. 127-142 


Aragón-Hernández € Bladé, Modelación numérica de flujo mixto en conductos cerrados con esquemas en volúmenes finitos 


Donde r, es un indicador del salto que sufre 
la solución a través de una onda en el contorno 
aguas arriba del elemento de volumen respecto 
del mismo salto en el otro contorno. Una posible 
expresión es el cociente entre contribuciones 
de la onda ¡ y el término que hace que el es- 
quema sea de segundo orden, ya que éstas son 
precisamente las que provocan las oscilaciones 
espurias, es: 


E (1 Ñ lv Dl... 


(27) 
EN (1 il. 


(; Jm > 


Por otra parte, para el tratamiento del térmi- 
no independiente H, es conveniente la descom- 
posición del mismo en función de los términos 
de fricción y de los términos de la pendiente del 
fondo, como: 


H=H!'+H?*= o (28) 
gL, +9 A(S, -S,) 
0 
1 = 

| 81,+34S, pd 

0 
de =gAS 30) 

f 


En la ecuación (10), H' es la expresión numé- 
rica del término independiente H y representa 
el término independiente integrado en todo el 
volumen finito; por lo tanto, considerando (28), 
se tiene: 


4,=H 2H (31) 
El término H?, que incluye los términos de 


fricción, se puede considerar en una discretiza- 
ción centrada simple como: 


H¡*=Ax,H?" (32) 
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Mientras que según Vázquez-Cendón (1999), 
H?, los términos de la pendiente del fondo 
deben discretizarse de acuerdo con el esquema 
numérico utilizado; para satisfacer la propiedad 
de conservación exacta puede ser dividido en 
dos contribuciones en las fronteras del volumen 
finito y descompuesto sobre los vectores propios 
de la matriz jacobiana J dela siguiente manera: 


He == a + Hao (33) 


Donde: 


aa = 50 +stno(»)fis(] (34) 


¡+1/2 


Los coeficientes o los incrementos de 
calado o carga de presión Ah, de las fuerzas 
de presión Al, y de la elevación del fondo del 
colector Az quedan definidos como: 


B, 58 Á(42+A1)+28A1, (35) 
B,=B, 

Ai =h,.-h, (36) 

Al, =1,, 1, (37) 

Az =Z,,-2, (38) 


Los esquemas numéricos como los anteriores 
son esquemas explícitos que tienen un coste 
computacional pequeño en cada paso de tiem- 
po, pero para ser estables necesitan trabajar con 
incrementos de tiempo también pequeños y de- 
ben cumplir con la condición de estabilidad de 
Courant-Friedrichs-Lewy (Chaudhry, 1993). La 
condición significa que ninguna onda presente 
en el dominio viaja más de una distancia Ax en 
un tiempo Af. Así, para el esquema numérico 
se tiene: 


(39) 
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Lo cual también se pueden escribir según el 
llamado número de Courant-Friedrichs-Lewy 
(CFL), como: 


Ju =c|At 
= LA <1l 


CEL 0 (40) 


Validación numérica 


Para validar el modelo numérico desarrollado, 
en esta sección se presentan un caso de referen- 
cia y dos ensayos de laboratorio existentes en 
la literatura técnica. En los casos de validación 
se presentan flujos que son muy comunes que 
ocurran en colectores pluviales y otros conduc- 
tos cerrados. Los tipos de flujo involucrados son 
flujo en lámina libre y flujo mixto. En el caso 
de flujo en lámina libre, puede existir flujo con 
régimen lento, rápido y transcrítico (de lento a 
rápido y de rápido a lento). 

Así, el flujo en régimen lento aparece en los 
tres casos de validación numérica, mientras que 
el flujo en régimen rápido se exhibe en el primer 
y tercer caso de validación numérica; por otro 
lado, el flujo en régimen transcrítico de lento a 
rápido ocurre en el primer caso de validación; 
mientras que, de rápido a lento a través de un 
resalto hidráulico, se manifiesta en el tercer en- 
sayo y, en este mismo, se presenta flujo mixto. 
Dichos casos se exponen a continuación. 


Flujo en régimen transcrítico 


Para demostrar la habilidad del modelo para re- 
producir régimen transcrítico, éste es aplicado a 
un caso de referencia sencillo; este caso sirve pa- 
ra validar el desempeño del esquema numérico 
al discretizar los términos convectivos, que por 
lo general causan las inestabilidades numéricas; 
es decir, no se consideran la influencia de la 
pendiente del fondo y la pendiente de fricción. 

Para ello, se toma en cuenta un canal pris- 
mático horizontal con una sobreelevación suave 
del fondo de forma parabólica. Utilizando las 
condiciones de contorno adecuadas se consigue 
que el flujo sobre el obstáculo sea en régimen 
lento, rápido o transcrítico (de rápido a lento 


o de lento a rápido). En esta oportunidad se 
considera flujo en régimen transcrítico de lento 
a rápido. El canal de sección rectangular tiene 
1 000 m de largo, 10 m de ancho y sin rugosidad. 
Las condiciones de contorno son las siguientes: 
aguas arriba, un caudal constante de 300 n?/s; 
aguas abajo se considera flujo supercrítico 
(cuando aguas abajo de la sobreelevación se 
presenta flujo subcrítico, se utiliza la condición 
de calado crítico; es decir, existe una caída libre). 
La condición inicial es un caudal nulo y la su- 
perficie libre del agua correspondiente al calado 
crítico impuesto aguas abajo. La cota del fondo 
del canal viene dada por: 


3-2 (x-500) 200m<x<800m 
3000 000 
z(x) = si (41) 
0 200m=x>800m 


El dominio de cálculo se discretizó en celdas 
de 4 m. De esta forma, la cota de la lámina de 
agua resultado de la aplicación del modelo se 
observa en la figura 2. El flujo crítico se presenta 
cuando el número de Froude es igual a uno; esto 
se produce en el punto localizado en x = 500 m. 
En este punto, el calado corresponde al calado 
crítico con un valor de 4.51 m. Además, en este 
mismo punto se produce el cambio de régimen, 
aguas arriba del mismo existe régimen lento y 
aguas abajo régimen rápido. 

También, en la figura 3 se pueden ver los va- 
lores de la velocidad a lo largo del canal; dichos 
valores son suaves y continuos. Además, y para 
mostrar la bondad del método implementado, se 
utilizó el método paso a paso (Chaudhry, 1993) 
para calcular la curva de remanso correspon- 
diente. Los resultados obtenidos con el modelo 
numérico y el calculado con el método paso a 
paso son muy similares, por lo que los obtenidos 
con este último no se graficaron debido a que las 
diferencias son imperceptibles, siendo siempre 
menores a 1 mm. Tanto los valores del calado 
como los de las velocidades, en conjunto con la 
comparación mencionada, dejan evidencia de la 
capacidad del esquema numérico para simular 
flujo en régimen permanente. 
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== Fondo del canal === + Lámina de agua 


Cota (m) 


0 100 200 300 400 500 600 700 800 900 1 000 
Longitud (m) 


Figura 2. Cota de la lámina de agua y fondo del canal. 
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'S Propagación hidráulica en un colector se realizaron en un colector circular con diá- 
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semejanza de Froude (White, 2004). Siguiendo 
a León, Ghidaoui, Schmidt y García (2006), 
al igual que a Sanders y Bradford (2011), los 
datos escalados de laboratorio son una tube- 
ría con una longitud de 304.8 m, diámetro de 
0.3048 m, pendiente de 0.001, coeficiente de 
fricción de Manning de 0.0115 s/m!” para 
el caudal base y 0.0116 s/m!”* para el caudal 
punta. La prueba se inicia en flujo permanen- 
te, con un caudal base de 0.004984 m?*/s y un 
calado normal en el extremo aguas abajo de 
0.0768 m; después se impone un hidrograma 
trapecial simétrico. El hidrograma es caracte- 
rizado por un caudal punta de 0.018689 m*/s 
durante 12 s, y 60 s en las ramas ascendente y 
descendente. Aguas abajo se consideró caída 
libre para forzar la existencia de flujo crítico. 
Los calados fueron medidos en dos puntos 
localizados en x = 8.66 m y x = 77.94 m desde 
el extremo aguas arriba. 

El modelo numérico de flujo en lámina libre 
se aplicó utilizando 348 volúmenes finitos. En 
la figura 4 se observan los calados medidos y 
calculados en x = 8.66 m y x = 77.94 m. 


La figura 4 muestra que el modelo predice 
adecuadamente el hidrograma en el punto 
localizado en x = 77.94 m, la rama ascendente y 
el calado en la punta del hidrograma del punto 
correspondiente a x = 8.66 m; sin embargo, en la 
rama descendente existen algunas diferencias. 
Las diferencias mostradas fueron reportadas 
también por otros autores (León et al., 2006; 
Sanders éz Bradford, 2011), quienes reprodu- 
jeron numéricamente este experimento. Según 
León et al. (2006), con base en otros autores, las 
diferencias son consecuencia de posibles proble- 
mas de escala y al hecho de que el experimento 
incluye caudales en la región de transición entre 
flujo laminar y flujo completamente turbulento. 
Lo anterior causa que el coeficiente de fricción 
de Mamning sea una función del número de 
Reynolds. 


Entrada en carga desde aguas abajo 
En la literatura técnica, el flujo mixto se ha 


estudiado de forma experimental y numérica, 
en especial con la entrada en carga desde aguas 


Calado (m) 


. Medido x = 8.66 m 
5 Medido x = 77.94 m 
Calculado x = 8.66 m 
= = - Calculado x = 77.94 m 


500 600 700 800 900 1000 
Tiempo (s) 


Figura 4. Calados medidos y calculados en x = 8.66 m y x =77.94 m. 
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abajo; pero el flujo mixto con entrada en carga 
desde aguas arriba y por ambos extremos de 
forma simultánea también es posible. En este 
trabajo se modela numéricamente el experi- 
mento de flujo mixto de Trajkovic et al. (1999), 
correspondiente a la entrada en carga desde 
aguas abajo. 

El dispositivo experimental utilizado con- 
siste de una tubería de metacrilato de sección 
circular de 10 m de longitud, un diámetro 
interno de 0.10 m y un coeficiente de rugosidad 
de Manmning de 0.008 s/m!'”. Se colocaron dos 
compuertas automáticas en los extremos aguas 
arriba y aguas abajo de la tubería, respectiva- 
mente. Con la operación de estas compuertas, 
las condiciones del flujo se pueden cambiar de 
forma rápida. Se colocaron tubos de ventilación 
en la clave de la tubería para evitar la posible 
interferencia de la fase de aire. Se llevaron a 
cabo diferentes tipos de ensayos, variando la 
pendiente de la tubería, condiciones iniciales y 
condiciones de contorno a través de la operación 
de las compuertas (Trajkovic et al., 1999). 

En los experimentos tipo A, los autores utili- 
zaron una pendiente de la tubería del 2.7%. Las 
condiciones del experimento fueron un caudal 
de entrada constante de 0.0013 m?/s, correspon- 
diente a una apertura de la compuerta en el ex- 
tremo aguas arriba de e, =0.014 m, y en el extre- 
mo aguas abajo, la compuerta totalmente abierta 
con un calado normal de 0.1 veces el diámetro. 
Debido a ello, el tipo de flujo es permanente y 
supercrítico. En un instante de tiempo dado del 
experimento, la compuerta del extremo aguas 
abajo fue cerrada rápidamente (cierre no instan- 
táneo), generando un transitorio en forma de 
frente de onda positivo moviéndose hacia aguas 
arriba. Después de 30 s del cierre, la compuerta 
fue reabierta de forma parcial, produciendo otro 
fenómeno transitorio. Se probaron diferentes 
valores para la reapertura de la compuerta; el 
correspondiente a una apertura de e, = 0.008 m 
es el considerado para la validación. Con esta 
apertura se registró una pequeña disminución 
de la carga de presión en el instante de tiempo 
30 s, pero debido a que el caudal de salida en 
la tubería es menor al de entrada, la carga de 


presión continuó incrementándose. Se midieron 
valores de la carga de presión en lo puntos P5 y 
P7, localizados en x =7.2 m y x= 9.2 m desde el 
extremo aguas arriba, respectivamente. 

Para la modelación numérica de este expe- 
rimento se utilizaron volúmenes finitos de 0.10 
m, una celeridad de la onda de 6 m/s, corres- 
pondiente a una ranura de Preissmann de 0.021 
veces el diámetro; el valor de la celeridad se 
utilizó como una variable de calibración, por lo 
que corresponde a la que tiene un error relativo 
menor. La carga de presión medida y calculada 
se muestra en la figura 5. 

En la figura 5 se puede observar que la carga 
de presión calculada coincide con la medida, 
incluso cuando se presenta la reapertura de 
la compuerta; el error relativo es de 3.7 y 3.1% 
para el punto P5 y P7, respectivamente. Por 
otra parte, la velocidad del frente de onda fue 
reproducida con buena aproximación (dife- 
rencia menor a 1 s en el punto P7 e inexistente 
en el punto P5); en cambio, la forma del frente 
de onda calculado es más pronunciada que el 
medido (diferencias máximas de 2.3 y 2.7 cm en 
los puntos P5 y P7, respectivamente). Algunos 
autores (León et al., 2009; Trajkovic et al., 1999) 
que han reproducido este experimento asumen 
que esta diferencia se debe a que en los cálculos 
el cierre de la compuerta aguas abajo se conside- 
ra de manera instantánea, mientras que en los 
ensayos hay un periodo de tiempo para realizar 
tal acción. En este trabajo se simuló el tiempo de 
cierre de la compuerta con distintas duraciones 
(0-2 s); los resultados no mejoraron la forma 
del frente de onda, por lo que dicha diferencia 
podría tener otro origen, probablemente rela- 
cionada con la discontinuidad del flujo. Para 
los esquemas numéricos empleados, el salto 
hidráulico es una discontinuidad en la solución 
y se modela como tal (las variables hidráulicas 
de la celda dependen de las celdas vecinas 
solamente), y la discontinuidad se produce en 
una distancia corta, mientras que la longitud del 
resalto hidráulico que se produce en el ensayo es 
mayor. Lo anterior podría mejorarse empleando 
volúmenes finitos de menor tamaño (acción no 
realizada). 
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4  Medido P7 +  Medido P5 
Calculado P7 


Calculado P5 


Clave del colector 


e 


Carga de presión (m) 
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Figura 5. Carga de presión medida y calculada en las secciones P5 y P7 para la reapertura de compuerta e, = 0.008 m. 


En la figura 6 se muestra la posición con 
respecto al tiempo y la velocidad del frente 
de onda a lo largo de la tubería. En ésta se 
observa que instantes posteriores a la apertura 
de la compuerta, la velocidad alcanza valores 
de hasta de 0.4 m/s; después, a medida que el 
resalto hidráulico se propaga hacia aguas arriba, 
el valor de la velocidad comienza a descender 
hasta llegar a un valor de 0.16 m/s al final de la 
simulación. Tanto la velocidad como la posición 
tienen un cambio en su comportamiento al mo- 
mento de la reapertura parcial de la compuerta, 
lo que provoca una disminución de la velocidad 
del frente de onda y, por tanto, un cambio de 
pendiente en la línea que representa la posición 
del mismo. 

Finalmente, en la figura 7 se muestran las 
condiciones del flujo antes del cierre de la com- 
puerta en el extremo aguas abajo del colector; el 
flujo es en lámina libre con régimen permanente 
y supercrítico. En el momento en que se reali- 
za el cierre rápido de la compuerta se genera 
un frente de onda que propicia un cambio de 
régimen (de rápido a lento), lo cual obliga a la 
formación de un resalto hidráulico que de inme- 
diato provoca la entrada en carga del colector y 
éste se propaga hacia aguas arriba. De la misma 
manera, en la reapertura parcial de la compuer- 
ta se aprecia la disminución de la velocidad del 


frente de onda, la cual se manifiesta como una 
disminución en la separación entre las líneas 
punteadas que representan la carga de presión 
a lo largo del colector en distintos instantes de 
tiempo, con una separación temporal constante 
(cada 10 s). 


Conclusiones 


En este trabajo se desarrolló un modelo nu- 
mérico para reproducir flujo en lámina libre 
y flujo en presión de forma simultánea (flujo 
mixto) en conductos cerrados. El modelo es 
capaz de reproducir flujo en lámina libre en 
régimen lento, rápido y transcrítico, con pre- 
sencia de discontinuidades sin necesidad de 
cálculos adicionales. La capacidad del modelo 
para simular flujo en régimen lento y rápido de 
forma automática es gracias a que el esquema 
numérico utilizado es un esquema conservativo 
y descentrado en volúmenes finitos. El flujo en 
presión se modela con la ayuda del método de 
la ranura de Preissmann. 

En el primer caso de la validación numérica, 
los resultados obtenidos permiten deducir que 
el modelo, en primer lugar, es capaz de repro- 
ducir flujo en régimen transcrítico de lento a 
rápido de forma precisa, tal como lo muestra 
la comparación con el método paso a paso. La 
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Figura 6. Posición y velocidad del frente de onda de presión. 
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Figura 7. Carga de presión en diferentes instantes de tiempo. 


localización y el valor del calado crítico son 
correctos, además de obtenerse diferencias en 
lámina de agua menores a 0.001 m, consideradas 
despreciables. En segundo lugar, el modelo es 
capaz de simular flujo permanente, que algunos 
modelos numéricos son incapaces de preservar. 
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En el segundo caso, el modelo se utilizó para 
la propagación hidráulica de un hidrograma en 
un colector circular en régimen lento. Se trata 
de una actividad muy común en la ingeniería 
práctica para determinar la atenuación y trasla- 
ción del hidrograma a su paso por un conducto. 
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Los resultados muestran que, de forma general, 
el modelo reproduce muy bien el experimento, 
salvo la rama descendente del punto de medida 
más próximo al extremo aguas arriba, x = 8.66 
m. Resultados similares fueron obtenidos por 
otros autores, indicando que las diferencias son 
debidas, entre otras causas, al problema del 
escalado de los datos. 

Finalmente, el tercer caso es muy especial, 
pues se presenta flujo en régimen supercrítico 
cuando el colector trabaja con la compuerta 
aguas abajo totalmente abierta; flujo en régimen 
transcrítico de rápido a lento a través de un re- 
salto hidráulico al cerrarse la compuerta aguas 
abajo, el cual se conserva durante el tiempo que 
dura el experimento; flujo en presión provocado 
también por el cierre de la compuerta aguas 
abajo, al alcanzar la lámina de agua la clave 
del colector. La variedad de flujos y fenómenos 
hidráulicos que se presentan en un espacio tan 
pequeño y un tiempo tan corto son un reto para 
los modelos numéricos; en este caso, el modelo 
numérico desarrollado fue capaz de reproducir 
el ensayo de forma más que aceptable y robusta. 
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